{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"../common/rfsoc_book_banner.jpg\" alt=\"University of Strathclyde\" align=\"left\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook Set B\n",
    "\n",
    "---\n",
    "\n",
    "## 03 - ADCs and DACs\n",
    "This notebook explores topics relating to the analogue-digital interface. These include aliasing and anti-aliasing filters for analogue to digital conversion and zero-order-hold techniques and reconstruction filters for digital to analogue conversion.\n",
    "\n",
    "## Table of Contents\n",
    "* [1. Introduction](#introduction)\n",
    "* [2. Analogue to Digital Conversion](#adc)\n",
    "    * [2.1. Aliasing](#aliasing)\n",
    "    * [2.2. Anti-Aliasing Filter](#anti-aliasing-filter)\n",
    "* [3. Digital to Analogue Conversion](#dac)\n",
    "    * [3.1. Zero-Order-Hold](#zero-order-hold)\n",
    "    * [3.2. Reconstruction Filter](#reconstruction-filter)\n",
    "* [4. Upper Nyquist Zones](#upper-nyquist-zones)\n",
    "* [5. Conclusion](#conclusion)\n",
    "\n",
    "## Revision\n",
    "* **v1.0** | 22/01/23 | *First Revision*\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction <a class=\"anchor\" id=\"introduction\"></a>\n",
    "Before we begin exploring the analogue-digital conversion process, we will load a few Python libraries below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import rfsoc_book.helper_functions as hf\n",
    "from scipy import signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Analogue to Digital Conversion <a class=\"anchor\" id=\"adc\"></a>\n",
    "We have already investigated two topics relating to analogue-digital conversion, which are sampling and quantisation. Aliasing is also of key importance when sampling an analogue signal, as higher frequency components (above the Nyquist sampling rate) can interfere with signals of interest.\n",
    "\n",
    "An anti-aliasing filter is typically used to suppress the effects of aliasing, as shown in Figure 1. Depending on the sampling frequency of the ADC, it may be possible to use a combination of oversampling and a cheap analogue filter to acquire the frequency band of interest.\n",
    "\n",
    "<figure>\n",
    "<img src=\"./images/analogue_anti_alias.png\" style=\"width: 75%;\"/>\n",
    "    <figcaption><b>Figure 1: Anti-aliasing filter placed before the Analogue to Digital Converter (ADC) to suppress the effects of aliasing. The top plot requires an expensive analogue filter to prevent aliasing. The bottom plot demonstrates the effects of oversampling the signal to use a cheaper analogue filter.</b></figcaption>\n",
    "</figure>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1. Aliasing <a class=\"anchor\" id=\"aliasing\"></a>\n",
    "We can easily demonstrate the effects of aliasing by sampling a signal that contains two frequency components. One frequency component will be within the Nyquist sampling rate and the other will be above the Nyquist sampling rate.\n",
    "\n",
    "Consider the cell below, which approximates a continuous signal that contains a 50Hz and 240Hz sine wave and some Gaussian noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Aa = 0.25   # Amplitude a\n",
    "Ab = 0.5    # Amplitude b\n",
    "fa = 50     # Desired frequency a\n",
    "fb = 240    # Desired frequency b\n",
    "t  = 0.2    # Time\n",
    "\n",
    "fs_cont = 2400    # Sampling rate of 'continuous' signal\n",
    "\n",
    "# Define 'continuous' signal for comparison\n",
    "x_cont = np.arange(0, t, 1/fs_cont)\n",
    "N_cont = np.size(x_cont)\n",
    "noise  = np.random.normal(0, 0.01, N_cont)\n",
    "y_cont = Aa * np.sin(2 * np.pi * fa * x_cont + 10) + Ab * np.sin(2 * np.pi * fb * x_cont) + noise\n",
    "\n",
    "# Plot the 'continuous' signal\n",
    "hf.plot_timeseries(\"Continuous Signal\",\n",
    "                   [x_cont], [y_cont], \n",
    "                   ['continuous'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the log-scale power spectra of the 'continuous' signal to easily inspect its frequency content. Notice that there are two lobes, one at 50Hz and the other at 240Hz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the log-scale power spectra of the 'continuous' signal\n",
    "Y_cont = np.fft.fft(y_cont)\n",
    "Y_cont_norm = np.abs(Y_cont)*2/N_cont\n",
    "Y_cont_log = 20*np.log10(Y_cont_norm)\n",
    "\n",
    "# Plot the log-scale power spectra\n",
    "fig = plt.figure(figsize=(10,4))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(np.arange(0, N_cont//2+1)*fs_cont/N_cont, Y_cont_log[0:N_cont//2+1])\n",
    "axes.set_title('Power Spectra of the Continuous Signal')\n",
    "axes.set_xlabel('Frequency (Hz)')\n",
    "axes.set_ylabel('Power Spectra (dB)')\n",
    "axes.grid(True, 'Major')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will sample the 'continuous' signal using a sampling frequency of 200Hz. According to Nyquist Sampling Theorem, all frequencies above half the sampling rate will alias. The sampled signal is plotted below in the time domain for inspection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sampled signal\n",
    "fs = 200\n",
    "x = np.arange(0, t, 1/fs)   # Discrete time i.e. sampled time period\n",
    "N = np.size(x)\n",
    "y = y_cont[0::fs_cont//fs]  # Get sampled signal\n",
    "\n",
    "# Plot the sampled signal\n",
    "hf.plot_timeseries(\"Sampled Signal\",\n",
    "                   [x_cont, x], [y_cont, y], \n",
    "                   ['continuous','discrete'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The log-scale power spectra of the sampled signal can be inspected below. Notice that the 50Hz frequency component is present in the output plot. However, there is also a 40Hz frequency component that was not in the original 'continuous' signal. This spurious frequency has emerged due to undersampling the 240Hz frequency component in the 'continuous' signal, which has caused it to alias into the sampled spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the sampled signal's log-scale power spectra\n",
    "Y = np.fft.fft(y)\n",
    "Y_norm = np.abs(Y)*2/N\n",
    "Y_log = 20*np.log10(Y_norm)\n",
    "\n",
    "# Plot the log-scale power spectra of the sampled signal\n",
    "fig = plt.figure(figsize=(10,4))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(np.arange(0, N//2+1)*fs/N, Y_log[0:N//2+1])\n",
    "axes.set_title('Power Spectra of the Sampled Signal')\n",
    "axes.set_xlabel('Frequency (Hz)')\n",
    "axes.set_ylabel('Power Spectra (dB)')\n",
    "axes.grid(True, 'Major')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The frequency of the aliased signal can be easily computed as shown in Figure 2. The sampled signal will fold across several Nyquist Zones before residing in Nyquist Zone 1. Once you have completed this notebook, you can try different frequencies as shown in the plot i.e. 180Hz and 310Hz.\n",
    "\n",
    "<figure>\n",
    "<img src=\"./images/nyquist_zones_examples.png\" style=\"width: 75%;\"/>\n",
    "    <figcaption><b>Figure 2: Examples of aliasing with reference to Nyquist Zones.</b></figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2. Anti-Aliasing Filter <a class=\"anchor\" id=\"anti-aliasing-filter\"></a>\n",
    "An anti-alias filter is used before the 'continuous' signal is sampled by the ADC. The purpose of the filter is to suppress frequencies above the Nyquist Sampling rate. We can generate a lowpass filter for this purpose below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A simple function to generate a lowpass filter\n",
    "def generate_lowpass(fs, cutoff, stop, numtaps):\n",
    "    taps = signal.remez(numtaps, [0, cutoff, stop, 0.5*fs], [1, 0], Hz=fs)\n",
    "    w, h = signal.freqz(taps, [1], worN=2000)\n",
    "    return w, h, taps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use the lowpass function to create a 151 tap filter and plot the response."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Obtain the anti-alias response and plot\n",
    "numtaps = 151\n",
    "w, h, coeffs = generate_lowpass(fs_cont, fs/2, 1.4*(fs/2), numtaps)\n",
    "hf.plot_response(fs_cont, w, h, title = \"Low-Pass Filter\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The lowpass filter is convolved with the 'continuous' signal in our simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter the 'continuous' signal\n",
    "y_cont_anti = signal.filtfilt(coeffs, 1, y_cont)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plot the filtered 'continuous' signal for inspection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the 'continuous' signal filtered by the anti-alias filter\n",
    "hf.plot_timeseries(\"Continuous Signal Filtered by Anti-Alias Filter\",\n",
    "                   [x_cont], [y_cont_anti],\n",
    "                   ['continuous'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that a signal in the time domain is a 50Hz sine wave. We can plot the log-scale power spectra of the signal to verify that the 240Hz frequency component has been suppressed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the log-scale power spectra of the anti-aliased signal\n",
    "Y_cont_anti = np.fft.fft(y_cont_anti)\n",
    "Y_cont_anti_norm = np.abs(Y_cont_anti)*2/N_cont\n",
    "Y_cont_anti_log = 20*np.log10(Y_cont_anti_norm)\n",
    "\n",
    "# Plot the log-scale power spectra of the anti-aliased signal\n",
    "fig = plt.figure(figsize=(10,4))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(np.arange(0, N_cont//2+1)*fs_cont/N_cont, Y_cont_anti_log[0:N_cont//2+1])\n",
    "axes.set_title('Power Spectra of the Continuous Signal Filtered by the Anti-Alias Filter')\n",
    "axes.set_xlabel('Frequency (Hz)')\n",
    "axes.set_ylabel('Power Spectra (dB)')\n",
    "axes.grid(True, 'Major')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It appears that the 240Hz sine wave is no longer present in the filtered 'continuous' signal. We can now sample the signal at 200Hz and plot the resulting log-scale power spectra below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample the continuous signal filtered by the anti-alias filter\n",
    "y_anti = y_cont_anti[0::fs_cont//fs]\n",
    "\n",
    "# Compute the log-scale power spectra of the sampled filtered signal\n",
    "Y_anti = np.fft.fft(y_anti)\n",
    "Y_anti_norm = np.abs(Y_anti)*2/N\n",
    "Y_anti_log = 20*np.log10(Y_anti_norm)\n",
    "\n",
    "# Plot the log-scale power spectra of the sampled filtered signal\n",
    "fig = plt.figure(figsize=(10,4))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(np.arange(0, N//2+1)*fs/N, Y_anti_log[0:N//2+1])\n",
    "axes.set_title('Power Spectra of the Sampled Wave Filtered by the Anti-Alias Filter')\n",
    "axes.set_xlabel('Frequency (Hz)')\n",
    "axes.set_ylabel('Power Spectra (dB)')\n",
    "axes.grid(True, 'Major')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 40Hz component seen previously is no longer present. The anti-aliasing filter has successfully successfully suppressed the 240Hz sine wave and prevented aliasing. Note that this example is for suppressing frequencies outside Nyquist Zone 1. When acquiring signals in higher order Nyquist Zones, bandpass filters should be used instead to suppress the effects of aliasing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Digital to Analogue Conversion <a class=\"anchor\" id=\"dac\"></a>\n",
    "Reconstruction is performed when converting a signal from the digital domain into analogue. A digital to analogue converter (DAC) will hold the voltage value for a period of $ t_s$ for each of the samples giving a \"steppy\" signal. An example of this process is illustrated in Figure 3.\n",
    "\n",
    "<figure>\n",
    "<img src=\"./images/DAC_reconstruction.png\" style=\"width: 75%;\"/>\n",
    "    <figcaption><b>Figure 3: DAC time domain signals and reconstruction filtering.</b></figcaption>\n",
    "</figure>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. Zero-Order-Hold\n",
    "To model the signal reconstruction performed by a DAC, we can apply a zero order hold (ZOH) to the digital signal. The function below acheives this, given an input signal and a sampling frequency. \n",
    "\n",
    "The zero order hold function below first makes the signal \"continuous\" by zero padding the array to a rate of 48000 Hz - the rate we use to approximate an analogue signal in this notebook. This padded signal is then convolved with the impulse response of the ZOH, which is a rectangular window of length, $ t_s $. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def zero_order_hold(signal, fs):\n",
    "    # New time axis for 'continuous' signal\n",
    "    x_zoh = np.arange(0, len(signal)/fs, 1/48000)\n",
    "    \n",
    "    # Zero pad the original signal to new rate\n",
    "    y_zoh = np.zeros(len(x_zoh))\n",
    "    y_zoh_indices = y_zoh[::int(np.round(48000/fs))]\n",
    "    y_zoh_indices[0:len(signal)] = signal\n",
    "    y_zoh[::int(np.round(48000/fs))] = y_zoh_indices\n",
    "    \n",
    "    # Create ZOH impulse response and convolve with input signal\n",
    "    h_zoh = np.ones(int(48000/fs))\n",
    "    y_zoh = np.convolve(h_zoh,y_zoh)\n",
    "    \n",
    "    # Chop of end to match time axis length\n",
    "    y_zoh = y_zoh[0:len(x_zoh)] \n",
    "    \n",
    "    return x_zoh, y_zoh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try changing the sampling rate to see the different outputs from the ZOH."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = 1       # Amplitude\n",
    "f = 100     # Desired frequency \n",
    "t = 0.2    # Time\n",
    "fs = 80    # Sampling rate\n",
    "\n",
    "# Define 'continuous' signal for comparison\n",
    "x_cont = np.arange(0, t, 1/48000)            \n",
    "y_cont = A * np.sin(2 * np.pi * f * x_cont + 10)  \n",
    "\n",
    "# Sampled signal\n",
    "x = np.arange(0, t, 1/fs)                # Discrete time i.e. sampled time period\n",
    "y = A * np.sin(2 * np.pi * f * x + 10)   # Formula for sine wave\n",
    "\n",
    "hf.plot_timeseries(\"Sampled Signal\",\n",
    "                   [x_cont, x], [y_cont, y], \n",
    "                   ['continuous','discrete'])\n",
    "\n",
    "# Perform ZOH on sampled signal\n",
    "[x_zoh, y_zoh] = zero_order_hold(y, fs)\n",
    "ax = hf.plot_timeseries(\"Zero Order Hold\", \n",
    "                        [x_zoh], [y_zoh], \n",
    "                        ['continuous'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. Reconstruction Filter <a class=\"anchor\" id=\"reconstruction-filter\"></a>\n",
    "\n",
    "Note that the output above is a little “steppy”, which is caused by the zero order hold (step reconstruction). This artefact can be removed with a reconstruction filter, which is implemented after the DAC using analogue circuitry. The reconstruction filter removes the baseband image and high frequency components present in the signal (in the form of the steps between the discrete levels).\n",
    "\n",
    "Ideally reconstruction filters have very sharp cut-off filters at $f_{s}/2$. Steeper roll-offs are more expensive, but clearly for many applications, good analogue filters are essential. In a DSP system, the precisely trimmed analogue filters could actually be more costly than the other DSP components, and their accuracy can be affected by temperature, ageing, etc. Therefore, it is desirable to do as much of the filtering as possible in the digital domain.\n",
    "\n",
    "To this end, we will first interpolate a sampled signal digitally, before passing it into our \"analogue\" reconstruction filter. Interpolation is simply upsampling a signal followed by a lowpass anti-imaging filter. Let us first design the interpolation filter below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = 3\n",
    "fs_cont = 48000\n",
    "Fs = fs*R\n",
    "w, h, coeffs = generate_lowpass(Fs, fs/2, 1.4*(fs/2), 15)\n",
    "hf.plot_response(Fs, w, h, title = \"Low-Pass Filter\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now apply the interpolation filter to the discrete signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_ov = np.arange(0, t, 1/(Fs))\n",
    "y_ov = np.zeros(len(y)*R)\n",
    "y_ov[::R] = y\n",
    "\n",
    "y_lp = signal.filtfilt(coeffs, 1, y_ov) * R\n",
    "hf.plot_timeseries(\"Interpolated Signal, R = 3\",\n",
    "                   [x_ov], [y_lp],\n",
    "                  ['discrete'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Zero Order Hold is then applied to the interpolated signal to model the DAC sampling process. The result is plotted below for inspection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform ZOH on sampled signal\n",
    "[x_zoh, y_zoh] = zero_order_hold(y_lp, Fs)\n",
    "hf.plot_timeseries(\"Zero-Order-Hold\",\n",
    "                   [x_zoh], [y_zoh],\n",
    "                   ['continuous'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A reconstruction filter with 400 filter taps is designed, which simulates an analogue reconstruction filter that is situated after the DAC. The filter response can be visualised below. The passband is very small, only allowing a few frequencies to remain and the rest are suppressed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w, h, coeffs = generate_lowpass(fs_cont, fs/2, Fs/2, 400)\n",
    "hf.plot_response(fs_cont, w, h, title = \"Low-Pass Filter\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of the reconstruction filter can be plotted alongside the step response for comparison. This plot is generated running the code cell below. Notice that the output of the filter is smooth and does not contain abrupt changes in amplitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [x_zoh, y_zoh] = zero_order_hold(y, fs)\n",
    "y_recon = signal.filtfilt(coeffs, 1, y_zoh)\n",
    "hf.plot_timeseries(\"Reconstructed Signal and Zero-Order-Hold Signal\",\n",
    "                   [x_zoh, x_zoh], [y_zoh, y_recon],\n",
    "                   ['continuous','continuous'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The smooth output response is more suitable for transmission as it does not contain spurious high frequency components that are present in the steppy output signal of the DAC."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Conclusion <a class=\"anchor\" id=\"conclusion\"></a>\n",
    "This notebook has simulated the analogue to digital conversion process and explored anti-aliasing filters. Additionally, the digital to analogue conversion process was investigated and reconstruction filters were discussed.\n",
    "\n",
    "In the next notebook, we review digital filters and design techniques."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "[⬅️ Previous Notebook](02_quantisation.ipynb) || [Next Notebook 🚀](04_digital_filter_design.ipynb)\n",
    "\n",
    "Copyright © 2023 Strathclyde Academic Media\n",
    "\n",
    "---\n",
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
